#ifndef _TemplateElement_
#define _TemplateElement_

#include "Dofs.h"
#include "GaussianPoint.h"
#include "Grid.h"
#include "Mesh.h"
#include <math.h>
using namespace std;
/**
 * @brief  建立一个参考单元，只考虑二维情况
 */
template <unsigned int DIM> class TemplateElement {
  public:
    /**
     * @brief 默认构造函数
     *
     */
    TemplateElement() = default;
    /**
     * @brief 默认析构函数
     *
     */
    virtual ~TemplateElement() = default;

    /**
     * @brief 从当前绑定的Grid中获取自由度信息
     *
     * @param i 指定序号；
     * @return Dofs<DIM>&　自由度的引用
     */
    virtual Dofs<DIM> &getDofsList(int i) = 0;
    /**
     * @brief 从当前绑定的Grid中获取第ｉ个自由度对应的全局编号；
     *
     * @param i 指定序号；
     * @return int 全局编号；
     */
    virtual int getGlobalIndex(int i) = 0;
    /**
     * @brief 返回本单元对应第i个高斯点坐标信息；
     *
     * @param i 指定序号
     * @return Point<DIM>& 高斯点
     */
    virtual Point<DIM> &GaussionPoint(int i) = 0;
    /**
     * @brief 返回本单元对应第i个高斯点权重信息；
     *
     * @param i 指定序号
     * @return double 高斯点对应权重的数值；
     */
    virtual double GaussionWeight(int i) = 0;
    /**
     * @brief 该类型单元对应的基函数；
     *
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i　基函数序数
     * @return double 该基函数的值；
     */
    virtual double phi(double xi, double eta, int i) = 0;
    /**
     * @brief 该类型单元对应的基函数对xi的偏导数；
     *
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i 序数
     * @return double 该偏导数的值
     */
    virtual double phi_xi(double xi, double eta, int i) = 0;
    /**
     * @brief 该类型单元对应的基函数对eta的偏导数；
     *
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i 序数
     * @return double 该偏导数的值
     */
    virtual double phi_eta(double xi, double eta, int i) = 0;
    /**
     * @brief
     * 任意单元到该模板单元仿射变换对应的雅克比矩阵行列式在(xi,eta)处的值；
     *
     * @param xi 坐标１
     * @param eta 坐标２
     * @return double 值
     */
    virtual double det_Jacobi(double xi, double eta) = 0;
    /**
     * @brief 基函数对x的偏导数；
     *
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i 序数
     * @return double 值
     */
    virtual double phi_x(double xi, double eta, int i) = 0;
    /**
     * @brief 基函数对y的偏导数；
     *
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i 序数
     * @return double 值
     */
    virtual double phi_y(double xi, double eta, int i) = 0;
    /**
     * @brief 计算并返回第ｉ个基函数在xi，eta点处的梯度向量；
     *
     * @param xi
     * @param eta
     * @param i
     * @return vector<double> 梯度向量；
     */
    virtual vector<double> gradient(double xi, double eta, int i) = 0;
    /**
     * @brief 标准单元中的点在真实单元的x坐标
     *
     * @param xi 坐标１
     * @param eta 坐标2
     * @return double 值
     */
    virtual double Global_x(double xi, double eta) = 0;
    /**
     * @brief 标准单元中的点在真实单元的y坐标
     *
     * @param xi 坐标１
     * @param eta 坐标2
     * @return double 值
     */
    virtual double Global_y(double xi, double eta) = 0;
    /**
     * @brief 该单元自由度个数；
     *
     * @return int
     */
    virtual int n_Dofs() = 0;
    /**
     * @brief 高斯点个数；
     *
     * @return int
     */
    virtual int n_GaussPnt() = 0;
    /**
     * @brief 计算并返回该单元的面积或体积；
     *
     * @return double
     */
    virtual double volume() = 0;
    /**
     * @brief 获取grid的自由度等信息
     * @details 将grid中的dofs赋给_dofslist
     */
    virtual void reinit(Grid<DIM> &_grid, Mesh<DIM> &_mesh) = 0;
    virtual void reinit(Grid<DIM> &_grid, Mesh<DIM> &_mesh, int index){};
    // virtual void setMesh(Mesh<DIM>& _mesh) = 0;
    // virtual vector<Dofs<2>> &getBoundaryDofsonGrid() = 0;

  protected:
    virtual double xi_x(double xi, double eta) = 0;
    virtual double xi_y(double xi, double eta) = 0;
    virtual double eta_x(double xi, double eta) = 0;
    virtual double eta_y(double xi, double eta) = 0;
};

/***************************************************************/
/********************** P1 ELEMENT *****************************/
/***************************************************************/

class Triangular_1_Element : public TemplateElement<2> {
    typedef vector<Dofs<2>> DofsList;
    typedef int DofsNum;

  public:
    Triangular_1_Element();
    ~Triangular_1_Element() = default;
    virtual Dofs<2> &getDofsList(int i);
    virtual int getGlobalIndex(int i);
    virtual Point<2> &GaussionPoint(int i) {
        return GaussInfo.GaussianPointInfo(i);
    };
    virtual double GaussionWeight(int i) {
        return GaussInfo.GaussianWeightInfo(i);
    };
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i);
    virtual double phi_eta(double xi, double eta, int i);
    virtual double det_Jacobi(double xi, double eta);
    virtual double xi_x(double xi, double eta);
    virtual double xi_y(double xi, double eta);
    virtual double eta_x(double xi, double eta);
    virtual double eta_y(double xi, double eta);
    virtual double phi_x(double xi, double eta, int i);
    virtual double phi_y(double xi, double eta, int i);
    virtual std::vector<double> gradient(double xi, double eta, int i);
    virtual double volume() { return det_Jacobi(0, 0) / 2.0; };
    virtual double Global_x(double xi, double eta);
    virtual double Global_y(double xi, double eta);
    virtual int n_Dofs() { return 3; };
    virtual int n_GaussPnt() { return GaussInfo.GaussianPointNumInfo(); };
    // TODO: new function
    virtual void reinit(Grid<2> &_grid, Mesh<2> &_mesh);
    virtual void reinit(Grid<2> &_grid, Mesh<2> &_mesh, int index);
    void reinit(Grid<2> &_grid, EasyMesh &easy_mesh, int index);

  private:
    DofsList dofslist;                         /*<dofs of grid*/
    Triangular_Element_GaussianInfo GaussInfo; /*<gauss information*/
};

Triangular_1_Element::Triangular_1_Element() {
    dofslist.resize(3);
    GaussInfo.LoadGaussianInfo(6);
}

Dofs<2> &Triangular_1_Element::getDofsList(int i) { return dofslist[i]; }

int Triangular_1_Element::getGlobalIndex(int i) {
    return dofslist[i].getGlobalIndex();
}

double Triangular_1_Element::phi(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return (1 - xi - eta);
    case 1:
        return xi;
    case 2:
        return eta;
    }
    return 0;
}

double Triangular_1_Element::phi_xi(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return -1;
    case 1:
        return 1;
    case 2:
        return 0;
    }
    return 0;
}

double Triangular_1_Element::phi_eta(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return -1;
    case 1:
        return 0;
    case 2:
        return 1;
    }
    return 0;
}

double Triangular_1_Element::det_Jacobi(double xi, double eta) {
    return (dofslist[1][0] - dofslist[0][0]) *
               (dofslist[2][1] - dofslist[0][1]) -
           (dofslist[1][1] - dofslist[0][1]) *
               (dofslist[2][0] - dofslist[0][0]);
}

double Triangular_1_Element::xi_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 3; i++)
        a = a + dofslist[i][1] * phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_1_Element::xi_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 3; i++)
        a = a + dofslist[i][0] * phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_1_Element::eta_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 3; i++)
        a = a + dofslist[i][1] * phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_1_Element::eta_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 3; i++)
        a = a + dofslist[i][0] * phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_1_Element::phi_x(double xi, double eta, int i) {
    return phi_xi(xi, eta, i) * xi_x(xi, eta) +
           phi_eta(xi, eta, i) * eta_x(xi, eta);
}

double Triangular_1_Element::phi_y(double xi, double eta, int i) {
    return phi_xi(xi, eta, i) * xi_y(xi, eta) +
           phi_eta(xi, eta, i) * eta_y(xi, eta);
}
std::vector<double> Triangular_1_Element::gradient(double xi, double eta,
                                                   int i) {
    std::vector<double> gradient({phi_x(xi, eta, i), phi_y(xi, eta, i)});
    return gradient;
}
double Triangular_1_Element::Global_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 3; i++) {
        a += dofslist[i][0] * phi(xi, eta, i);
    }
    return a;
}

double Triangular_1_Element::Global_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 3; i++) {
        a += dofslist[i][1] * phi(xi, eta, i);
    }
    return a;
}

void Triangular_1_Element::reinit(Grid<2> &_grid, Mesh<2> &_mesh) {
    // TODO:if structed mesh
}
void Triangular_1_Element::reinit(Grid<2> &_grid, Mesh<2> &_mesh, int index) {
    // TODO: if unstructed mesh, you must use this function,not prevous one.
    dofslist.clear();
    int num = this->n_Dofs();
    vector<Point<2>> p = _grid.getVertices();
    vector<pair<int, int>> dofsIndex = _mesh.getIndexofGrid(index);
    for (int k = 0; k < num; k++) {
        Dofs<2> q = {p[k], dofsIndex[k].first};
        dofslist.push_back(q);
    }
}

void Triangular_1_Element::reinit(Grid<2> &_grid, EasyMesh &easy_mesh,
                                  int index) {
    dofslist.clear();
    int num = this->n_Dofs();
    vector<Point<2>> p = _grid.getVertices();
    vector<pair<int, int>> dofsIndex = easy_mesh.getIndexofGrid(index);
    for (int k = 0; k < num; k++) {
        Dofs<2> q = {p[k], dofsIndex[k].first};
        dofslist.push_back(q);
    }
};

/***************************************************************/
/********************** P2 ELEMENT *****************************/
/***************************************************************/

class Triangular_2_Element : public TemplateElement<2> {
  public:
    Triangular_2_Element();
    ~Triangular_2_Element();
    virtual Dofs<2> &getDofsList(int i);
    virtual int getGlobalIndex(int i);
    virtual Point<2> &GaussionPoint(int i) {
        return GaussInfo.GaussianPointInfo(i);
    };
    virtual double GaussionWeight(int i) {
        return GaussInfo.GaussianWeightInfo(i);
    };
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i);
    virtual double phi_eta(double xi, double eta, int i);
    virtual double det_Jacobi(double xi, double eta);
    virtual double xi_x(double xi, double eta);
    virtual double xi_y(double xi, double eta);
    virtual double eta_x(double xi, double eta);
    virtual double eta_y(double xi, double eta);
    virtual double phi_x(double xi, double eta, int i);
    virtual double phi_y(double xi, double eta, int i);
    virtual std::vector<double> gradient(double xi, double eta, int i);
    virtual double volume() { return det_Jacobi(0, 0); };
    virtual double Global_x(double xi, double eta);
    virtual double Global_y(double xi, double eta);
    virtual int n_Dofs() { return 6; };
    virtual int n_GaussPnt() { return GaussInfo.GaussianPointNumInfo(); };
    // TODO: new function
    virtual void reinit(Grid<2> &_grid, Mesh<2> &_mesh);
    virtual void reinit(Grid<2> &_grid, Mesh<2> &_mesh, int index);
    void reinit(Grid<2> &_grid, EasyMesh &easy_mesh, int index);

  private:
    std::vector<Dofs<2>> dofslist;
    Triangular_Element_GaussianInfo GaussInfo;
};

Triangular_2_Element::Triangular_2_Element() {
    dofslist.resize(6);
    GaussInfo.LoadGaussianInfo(8);
}

Triangular_2_Element::~Triangular_2_Element(){};

Dofs<2> &Triangular_2_Element::getDofsList(int i) { return dofslist[i]; }
int Triangular_2_Element::getGlobalIndex(int i) {
    return dofslist[i].getGlobalIndex();
}
double Triangular_2_Element::phi(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return (1 - xi - eta) * (2 * (1 - xi - eta) - 1);
    case 1:
        return xi * (2 * xi - 1);
    case 2:
        return eta * (2 * eta - 1);
    case 5:
        return 4 * xi * (1 - xi - eta);
    case 3:
        return 4 * xi * eta;
    case 4:
        return 4 * eta * (1 - xi - eta);
    }
    return 0;
}

double Triangular_2_Element::phi_xi(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return 1 - 4 * (1 - xi - eta);
    case 1:
        return 4 * xi - 1;
    case 2:
        return 0;
    case 5:
        return 4 - 8 * xi - 4 * eta;
    case 3:
        return 4 * eta;
    case 4:
        return -4 * eta;
    }
    return 0;
}

double Triangular_2_Element::phi_eta(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return 1 - 4 * (1 - xi - eta);
    case 1:
        return 0;
    case 2:
        return 4 * eta - 1;
    case 5:
        return -4 * xi;
    case 3:
        return 4 * xi;
    case 4:
        return 4 - 4 * xi - 8 * eta;
    }
    return 0;
}

double Triangular_2_Element::det_Jacobi(double xi, double eta) {
    double J11 = 0, J12 = 0, J21 = 0, J22 = 0;
    for (int i = 0; i < n_Dofs(); i++) {
        J11 = J11 + dofslist[i][0] * phi_xi(xi, eta, i);
        J12 = J12 + dofslist[i][1] * phi_xi(xi, eta, i);
        J21 = J21 + dofslist[i][0] * phi_eta(xi, eta, i);
        J22 = J22 + dofslist[i][1] * phi_eta(xi, eta, i);
    }
    return (J11 * J22 - J12 * J21);
}

double Triangular_2_Element::xi_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++)
        a = a + dofslist[i][1] * phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_2_Element::xi_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++)
        a = a + dofslist[i][0] * phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_2_Element::eta_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++)
        a = a + dofslist[i][1] * phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_2_Element::eta_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++)
        a = a + dofslist[i][0] * phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_2_Element::phi_x(double xi, double eta, int i) {
    return phi_xi(xi, eta, i) * xi_x(xi, eta) +
           phi_eta(xi, eta, i) * eta_x(xi, eta);
}

double Triangular_2_Element::phi_y(double xi, double eta, int i) {
    return phi_xi(xi, eta, i) * xi_y(xi, eta) +
           phi_eta(xi, eta, i) * eta_y(xi, eta);
}
std::vector<double> Triangular_2_Element::gradient(double xi, double eta,
                                                   int i) {
    std::vector<double> gradient({phi_x(xi, eta, i), phi_y(xi, eta, i)});
    return gradient;
}
double Triangular_2_Element::Global_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++) {
        a += dofslist[i][0] * phi(xi, eta, i);
    }
    return a;
}

double Triangular_2_Element::Global_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++) {
        a += dofslist[i][1] * phi(xi, eta, i);
    }
    return a;
}
void Triangular_2_Element::reinit(Grid<2> &_grid, Mesh<2> &_mesh) {
    // TODO:if structed mesh
}
void Triangular_2_Element::reinit(Grid<2> &_grid, Mesh<2> &_mesh, int index) {
    // TODO: if unstructed mesh, you must use this function,not prevous one.
    dofslist.clear();
    int num = this->n_Dofs();
    vector<Point<2>> p = _grid.getVertices();
    vector<Point<2>> v = p;
    v.resize(6);
    // for (int j = 3; j < 6; j++){
    //    p[j] = (p[j-3] + p[(j-2)%3])/2.0; //computer midpoint as dofs
    // }
    v[3] = (v[1] + v[2]) / 2.0;
    v[4] = (v[0] + v[2]) / 2.0;
    v[5] = (v[0] + v[1]) / 2.0;
    vector<pair<int, int>> dofsIndex = _mesh.getIndexofGrid(index);
    for (int k = 0; k < num; k++) {
        Dofs<2> q = {p[k], dofsIndex[k].first}; // just add the vertices
        dofslist.push_back(q);
    }
}

void Triangular_2_Element::reinit(Grid<2> &_grid, EasyMesh &easy_mesh,
                                  int index) {
    dofslist.clear();
    int num = this->n_Dofs();
    // vector<Point<2>> p = _grid.getVertices();
    vector<pair<int, int>> dofsIndex = easy_mesh.getIndexofGrid(index);
    for (int k = 0; k < num; k++) {
        int index = dofsIndex[k].first;
        Dofs<2> dof = easy_mesh.getDofs(index);
        dofslist.push_back(dof);
    }
};
/***************************************************************/
/********************** P3 ELEMENT *****************************/
/***************************************************************/
class Triangular_3_Element : public TemplateElement<2> {
  public:
    Triangular_3_Element();
    ~Triangular_3_Element();
    // virtual void SetDofsList(std::vector<Dofs<2> >);
    virtual Dofs<2> &getDofsList(int i);
    // virtual int NdIdx(int i) ;
    virtual int getGlobalIndex(int i);
    virtual Point<2> &GaussionPoint(int i) {
        return GaussInfo.GaussianPointInfo(i);
    };
    virtual double GaussionWeight(int i) {
        return GaussInfo.GaussianWeightInfo(i);
    };
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i);
    virtual double phi_eta(double xi, double eta, int i);
    virtual double det_Jacobi(double xi, double eta);
    virtual double xi_x(double xi, double eta);
    virtual double xi_y(double xi, double eta);
    virtual double eta_x(double xi, double eta);
    virtual double eta_y(double xi, double eta);
    virtual double phi_x(double xi, double eta, int i);
    virtual double phi_y(double xi, double eta, int i);
    virtual std::vector<double> gradient(double xi, double eta, int i);
    virtual double volume() { return det_Jacobi(0, 0); };
    virtual double Global_x(double xi, double eta);
    virtual double Global_y(double xi, double eta);
    virtual int n_Dofs() { return 10; };
    virtual int n_GaussPnt() { return GaussInfo.GaussianPointNumInfo(); };
    // TODO: new function
    virtual void reinit(Grid<2> &_grid, Mesh<2> &_mesh);
    virtual void reinit(Grid<2> &_grid, Mesh<2> &_mesh, int index);
    void reinit(Grid<2> &_grid, EasyMesh &easy_mesh, int index);

  private:
    std::vector<Dofs<2>> dofslist;
    Triangular_Element_GaussianInfo GaussInfo;
};

Triangular_3_Element::Triangular_3_Element() {
    dofslist.resize(10);
    GaussInfo.LoadGaussianInfo(6);
}

Triangular_3_Element::~Triangular_3_Element(){};
/* delete
   void Triangular_3_Element::SetDofsList(std::vector<Dofs<2> > _list)
   {
   dofslist = _list;
   }
   */
Dofs<2> &Triangular_3_Element::getDofsList(int i) { return dofslist[i]; }
int Triangular_3_Element::getGlobalIndex(int i) {
    return dofslist[i].getGlobalIndex();
}

double Triangular_3_Element::phi(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return 4.5 * (1 - xi - eta) * (2.0 / 3 - xi - eta) *
               (1.0 / 3 - xi - eta);
    case 1:
        return 4.5 * xi * (xi - 1.0 / 3) * (xi - 2.0 / 3);
    case 2:
        return 4.5 * eta * (eta - 1.0 / 3) * (eta - 2.0 / 3);
    case 7:
        return 13.5 * xi * (1 - xi - eta) * (2.0 / 3 - xi - eta);
    case 8:
        return 13.5 * xi * (1 - xi - eta) * (xi - 1.0 / 3);
    case 3:
        return 13.5 * xi * eta * (xi - 1.0 / 3);
    case 4:
        return 13.5 * xi * eta * (eta - 1.0 / 3);
    case 5:
        return 13.5 * (1 - xi - eta) * eta * (eta - 1.0 / 3);
    case 6:
        return 13.5 * (1 - xi - eta) * (2.0 / 3 - xi - eta) * eta;
    case 9:
        return 27 * xi * eta * (1 - xi - eta);
    }
    return 0;
}

double Triangular_3_Element::phi_xi(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return -4.5 * (2.0 / 3 - xi - eta) * (1.0 / 3 - xi - eta) -
               4.5 * (1 - xi - eta) * (1.0 / 3 - xi - eta) -
               4.5 * (1 - xi - eta) * (2.0 / 3 - xi - eta);
    case 1:
        return 4.5 * (xi - 1.0 / 3) * (xi - 2.0 / 3) +
               4.5 * xi * (xi - 2.0 / 3) + 4.5 * xi * (xi - 1.0 / 3);
    case 2:
        return 0;
    case 7:
        return 13.5 * (xi + eta - 1) * (xi + eta - 2.0 / 3) +
               13.5 * xi * (xi + eta - 2.0 / 3) + 13.5 * xi * (xi + eta - 1);
    case 8:
        return 13.5 * (1 - xi - eta) * (xi - 1.0 / 3) -
               13.5 * xi * (xi - 1.0 / 3) + 13.5 * xi * (1 - xi - eta);
    case 3:
        return 13.5 * eta * (xi - 1.0 / 3) + 13.5 * xi * eta;
    case 4:
        return 13.5 * eta * (eta - 1.0 / 3);
    case 5:
        return -13.5 * eta * (eta - 1.0 / 3);
    case 6:
        return -13.5 * (2.0 / 3 - xi - eta) * eta - 13.5 * (1 - xi - eta) * eta;
    case 9:
        return 27 * eta * (1 - xi - eta) - 27 * xi * eta;
    }
    return 0;
}

double Triangular_3_Element::phi_eta(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return -4.5 * (2.0 / 3 - xi - eta) * (1.0 / 3 - xi - eta) -
               4.5 * (1 - xi - eta) * (1.0 / 3 - xi - eta) -
               4.5 * (1 - xi - eta) * (2.0 / 3 - xi - eta);
    case 1:
        return 0;
    case 2:
        return 4.5 * (eta - 1.0 / 3) * (eta - 2.0 / 3) +
               4.5 * eta * (eta - 2.0 / 3) + 4.5 * eta * (eta - 1.0 / 3);
    case 7:
        return 13.5 * xi * (xi + eta - 2.0 / 3) + 13.5 * xi * (xi + eta - 1);
    case 8:
        return -13.5 * xi * (xi - 1.0 / 3);
    case 3:
        return 13.5 * xi * (xi - 1.0 / 3);
    case 4:
        return 13.5 * xi * (eta - 1.0 / 3) + 13.5 * xi * eta;
    case 5:
        return -13.5 * eta * (eta - 1.0 / 3) +
               13.5 * (1 - xi - eta) * (eta - 1.0 / 3) +
               13.5 * (1 - xi - eta) * eta;
    case 6:
        return -13.5 * (2.0 / 3 - xi - eta) * eta -
               13.5 * (1 - xi - eta) * eta +
               13.5 * (1 - xi - eta) * (2.0 / 3 - xi - eta);
    case 9:
        return 27 * xi * (1 - xi - eta) - 27 * xi * eta;
    }
    return 0;
}

double Triangular_3_Element::det_Jacobi(double xi, double eta) {
    double J11 = 0, J12 = 0, J21 = 0, J22 = 0;
    for (int i = 0; i < n_Dofs(); i++) {
        J11 = J11 + dofslist[i][0] * phi_xi(xi, eta, i);
        J12 = J12 + dofslist[i][1] * phi_xi(xi, eta, i);
        J21 = J21 + dofslist[i][0] * phi_eta(xi, eta, i);
        J22 = J22 + dofslist[i][1] * phi_eta(xi, eta, i);
    }
    return (J11 * J22 - J12 * J21);
}

double Triangular_3_Element::xi_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++)
        a = a + dofslist[i][1] * phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_3_Element::xi_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++)
        a = a + dofslist[i][0] * phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_3_Element::eta_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++)
        a = a + dofslist[i][1] * phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_3_Element::eta_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++)
        a = a + dofslist[i][0] * phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_3_Element::phi_x(double xi, double eta, int i) {
    return phi_xi(xi, eta, i) * xi_x(xi, eta) +
           phi_eta(xi, eta, i) * eta_x(xi, eta);
}

double Triangular_3_Element::phi_y(double xi, double eta, int i) {
    return phi_xi(xi, eta, i) * xi_y(xi, eta) +
           phi_eta(xi, eta, i) * eta_y(xi, eta);
}
std::vector<double> Triangular_3_Element::gradient(double xi, double eta,
                                                   int i) {
    std::vector<double> gradient({phi_x(xi, eta, i), phi_y(xi, eta, i)});
    return gradient;
}
double Triangular_3_Element::Global_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++) {
        a += dofslist[i][0] * phi(xi, eta, i);
    }
    return a;
}

double Triangular_3_Element::Global_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++) {
        a += dofslist[i][1] * phi(xi, eta, i);
    }
    return a;
}
/* delete
   int Triangular_3_Element::NdIdx(int i)
   {

   return dofslist[i - 1].GlobalIndex();
   }
   */
void Triangular_3_Element::reinit(Grid<2> &_grid, Mesh<2> &_mesh) {
    // TODO:if structed mesh
}
void Triangular_3_Element::reinit(Grid<2> &_grid, Mesh<2> &_mesh, int index) {
    // TODO: if unstructed mesh, you must use this function,not prevous one.
    dofslist.clear();
    int num = this->n_Dofs();
    vector<Point<2>> p = _grid.getVertices();
    vector<Point<2>> v = p; // using v to store point
    v.resize(10);
    // for (int j = 1 ;j <= 3 ; j++){
    //   v[2*j+1] = (p[j-1] * 2.0 + p[j%3])/3.0; //trisection
    //   v[2*j+2] = (p[j-1] + p[j%3] * 2.0)/3.0; //trisection
    // }
    //   v[9] = (p[0] + p[1] + p[2])/3.0;
    v[7] = (v[0] * 2.0 + v[1]) / 3.0;
    v[8] = (v[0] + v[1] * 2.0) / 3.0;
    v[3] = (v[1] * 2.0 + v[2]) / 3.0;
    v[4] = (v[1] + v[2] * 2.0) / 3.0;
    v[5] = (v[2] * 2.0 + v[0]) / 3.0;
    v[6] = (v[2] + v[0] * 2.0) / 3.0;
    v[9] = (v[0] + v[1] + v[2]) / 3.0;
    vector<pair<int, int>> dofsIndex = _mesh.getIndexofGrid(index);
    for (int k = 0; k < num; k++) {
        Dofs<2> q = {v[k], dofsIndex[k].first};
        dofslist.push_back(q);
    }
}

void Triangular_3_Element::reinit(Grid<2> &_grid, EasyMesh &easy_mesh,
                                  int index) {
    dofslist.clear();
    int num = this->n_Dofs();
    vector<pair<int, int>> dofsIndex = easy_mesh.getIndexofGrid(index);
    for (int k = 0; k < num; k++) {
        int index = dofsIndex[k].first;
        Dofs<2> dof = easy_mesh.getDofs(index);
        dofslist.push_back(dof);
    }
};
/***************************************************************/
/********************** P4 ELEMENT *****************************/
/***************************************************************/
class Triangular_4_Element : public TemplateElement<2> {
  public:
    Triangular_4_Element();
    ~Triangular_4_Element();
    // virtual void SetDofsList(std::vector<Dofs<2> >);
    virtual Dofs<2> &getDofsList(int i);
    // virtual int NdIdx(int i) ;
    virtual int getGlobalIndex(int i);
    virtual Point<2> &GaussionPoint(int i) {
        return GaussInfo.GaussianPointInfo(i);
    };
    virtual double GaussionWeight(int i) {
        return GaussInfo.GaussianWeightInfo(i);
    };
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i);
    virtual double phi_eta(double xi, double eta, int i);
    virtual double det_Jacobi(double xi, double eta);
    virtual double xi_x(double xi, double eta);
    virtual double xi_y(double xi, double eta);
    virtual double eta_x(double xi, double eta);
    virtual double eta_y(double xi, double eta);
    virtual double phi_x(double xi, double eta, int i);
    virtual double phi_y(double xi, double eta, int i);
    virtual std::vector<double> gradient(double xi, double eta, int i);
    virtual double volume() { return det_Jacobi(0, 0); };
    virtual double Global_x(double xi, double eta);
    virtual double Global_y(double xi, double eta);
    virtual int n_Dofs() { return 15; };
    virtual int n_GaussPnt() { return GaussInfo.GaussianPointNumInfo(); };
    // TODO: new function
    virtual void reinit(Grid<2> &_grid, Mesh<2> &_mesh);
    virtual void reinit(Grid<2> &_grid, Mesh<2> &_mesh, int index);
    void reinit(Grid<2> &_grid, EasyMesh &easy_mesh, int index);

  private:
    std::vector<Dofs<2>> dofslist;
    Triangular_Element_GaussianInfo GaussInfo;
};

Triangular_4_Element::Triangular_4_Element() {
    dofslist.resize(15);
    GaussInfo.LoadGaussianInfo(8); //查表确定用哪个精度
}

Triangular_4_Element::~Triangular_4_Element(){};
/* delete
   void Triangular_4_Element::SetDofsList(std::vector<Dofs<2> > _list)
   {
   dofslist = _list;
   }
   */
Dofs<2> &Triangular_4_Element::getDofsList(int i) { return dofslist[i]; }
int Triangular_4_Element::getGlobalIndex(int i) {
    return dofslist[i].getGlobalIndex();
}

double Triangular_4_Element::phi(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return 32.0 / 3 * (0.25 - xi - eta) * (0.5 - xi - eta) *
               (0.75 - xi - eta) * (1 - xi - eta);
    case 1:
        return 32.0 / 3 * xi * (xi - 0.25) * (xi - 0.5) * (xi - 0.75);
    case 2:
        return 32.0 / 3 * eta * (eta - 0.25) * (eta - 0.5) * (eta - 0.75);
    case 9:
        return 128.0 / 3 * xi * (0.5 - xi - eta) * (0.75 - xi - eta) *
               (1 - xi - eta);
    case 10:
        return 64 * xi * (xi - 0.25) * (0.75 - xi - eta) * (1 - xi - eta);
    case 11:
        return 128.0 / 3 * xi * (xi - 0.25) * (xi - 0.5) * (1 - xi - eta);
    case 3:
        return 128.0 / 3 * xi * (xi - 0.25) * (xi - 0.5) * eta;
    case 4:
        return 64 * xi * (xi - 0.25) * eta * (eta - 0.25);
    case 5:
        return 128.0 / 3 * xi * eta * (eta - 0.25) * (eta - 0.5);
    case 6:
        return 128.0 / 3 * eta * (eta - 0.25) * (eta - 0.5) * (1 - xi - eta);
    case 7:
        return 64 * eta * (eta - 0.25) * (0.75 - xi - eta) * (1 - xi - eta);
    case 8:
        return 128.0 / 3 * eta * (0.5 - xi - eta) * (0.75 - xi - eta) *
               (1 - xi - eta);
    case 14:
        return 128 * xi * eta * (0.75 - xi - eta) * (1 - xi - eta);
    case 12:
        return 128 * xi * (xi - 0.25) * eta * (1 - xi - eta);
    case 13:
        return 128 * xi * eta * (eta - 0.25) * (1 - xi - eta);
    }
    return 0;
}

double Triangular_4_Element::phi_xi(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return -32.0 / 3 * (0.5 - xi - eta) * (0.75 - xi - eta) *
                   (1 - xi - eta) -
               32.0 / 3 * (0.25 - xi - eta) * (0.75 - xi - eta) *
                   (1 - xi - eta) -
               32.0 / 3 * (0.25 - xi - eta) * (0.5 - xi - eta) *
                   (1 - xi - eta) -
               32.0 / 3 * (0.25 - xi - eta) * (0.5 - xi - eta) *
                   (0.75 - xi - eta);
    case 1:
        return 32.0 / 3 * (xi - 0.25) * (xi - 0.5) * (xi - 0.75) +
               32.0 / 3 * xi * (xi - 0.5) * (xi - 0.75) +
               32.0 / 3 * xi * (xi - 0.25) * (xi - 0.75) +
               32.0 / 3 * xi * (xi - 0.25) * (xi - 0.5);
    case 2:
        return 0;
    case 9:
        return 128.0 / 3 * (0.5 - xi - eta) * (0.75 - xi - eta) *
                   (1 - xi - eta) -
               128.0 / 3 * xi * (0.75 - xi - eta) * (1 - xi - eta) -
               128.0 / 3 * xi * (0.5 - xi - eta) * (1 - xi - eta) -
               128.0 / 3 * xi * (0.5 - xi - eta) * (0.75 - xi - eta);
    case 10:
        return 64 * (xi - 0.25) * (0.75 - xi - eta) * (1 - xi - eta) +
               64 * xi * (0.75 - xi - eta) * (1 - xi - eta) -
               64 * xi * (xi - 0.25) * (1 - xi - eta) -
               64 * xi * (xi - 0.25) * (0.75 - xi - eta);
    case 11:
        return 128.0 / 3 * (xi - 0.25) * (xi - 0.5) * (1 - xi - eta) +
               128.0 / 3 * xi * (xi - 0.5) * (1 - xi - eta) +
               128.0 / 3 * xi * (xi - 0.25) * (1 - xi - eta) -
               128.0 / 3 * xi * (xi - 0.25) * (xi - 0.5);
    case 3:
        return 128.0 / 3 * (xi - 0.25) * (xi - 0.5) * eta +
               128.0 / 3 * xi * (xi - 0.5) * eta +
               128.0 / 3 * xi * (xi - 0.25) * eta;
    case 4:
        return 64 * (xi - 0.25) * eta * (eta - 0.25) +
               64 * xi * eta * (eta - 0.25);
    case 5:
        return 128.0 / 3 * eta * (eta - 0.25) * (eta - 0.5);
    case 6:
        return -128.0 / 3 * eta * (eta - 0.25) * (eta - 0.5);
    case 7:
        return -64 * eta * (eta - 0.25) * (1 - xi - eta) -
               64 * eta * (eta - 0.25) * (0.75 - xi - eta);
    case 8:
        return -128.0 / 3 * eta * (0.75 - xi - eta) * (1 - xi - eta) -
               128.0 / 3 * eta * (0.5 - xi - eta) * (1 - xi - eta) -
               128.0 / 3 * eta * (0.5 - xi - eta) * (0.75 - xi - eta);
    case 14:
        return 128 * eta * (0.75 - xi - eta) * (1 - xi - eta) -
               128 * xi * eta * (1 - xi - eta) -
               128 * xi * eta * (0.75 - xi - eta);
    case 12:
        return 128 * eta * (xi - 0.25) * (1 - xi - eta) +
               128 * xi * eta * (1 - xi - eta) - 128 * xi * eta * (xi - 0.25);
    case 13:
        return 128 * eta * (eta - 0.25) * (1 - xi - eta) -
               128 * xi * eta * (eta - 0.25);
    }
    return 0;
}

double Triangular_4_Element::phi_eta(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return -32.0 / 3 * (0.5 - xi - eta) * (0.75 - xi - eta) *
                   (1 - xi - eta) -
               32.0 / 3 * (0.25 - xi - eta) * (0.75 - xi - eta) *
                   (1 - xi - eta) -
               32.0 / 3 * (0.25 - xi - eta) * (0.5 - xi - eta) *
                   (1 - xi - eta) -
               32.0 / 3 * (0.25 - xi - eta) * (0.5 - xi - eta) *
                   (0.75 - xi - eta);
    case 1:
        return 0;
    case 2:
        return 32.0 / 3 * (eta - 0.25) * (eta - 0.5) * (eta - 0.75) +
               32.0 / 3 * eta * (eta - 0.5) * (eta - 0.75) +
               32.0 / 3 * eta * (eta - 0.25) * (eta - 0.75) +
               32.0 / 3 * eta * (eta - 0.25) * (eta - 0.5);
    case 9:
        return -128.0 / 3 * xi * (0.75 - xi - eta) * (1 - xi - eta) -
               128.0 / 3 * xi * (0.5 - xi - eta) * (1 - xi - eta) -
               128.0 / 3 * xi * (0.5 - xi - eta) * (0.75 - xi - eta);
    case 10:
        return -64 * xi * (xi - 0.25) * (1 - xi - eta) -
               64 * xi * (xi - 0.25) * (0.75 - xi - eta);
    case 11:
        return -128.0 / 3 * xi * (xi - 0.25) * (xi - 0.5);
    case 3:
        return 128.0 / 3 * xi * (xi - 0.25) * (xi - 0.5);
    case 4:
        return 64 * xi * (xi - 0.25) * (eta - 0.25) +
               64 * xi * (xi - 0.25) * eta;
    case 5:
        return 128.0 / 3 * xi * (eta - 0.25) * (eta - 0.5) +
               128.0 / 3 * xi * eta * (eta - 0.5) +
               128.0 / 3 * xi * eta * (eta - 0.25);
    case 6:
        return 128.0 / 3 * (eta - 0.25) * (eta - 0.5) * (1 - xi - eta) +
               128.0 / 3 * eta * (eta - 0.5) * (1 - xi - eta) +
               128.0 / 3 * eta * (eta - 0.25) * (1 - xi - eta) -
               128.0 / 3 * eta * (eta - 0.25) * (eta - 0.5);
    case 7:
        return 64 * (eta - 0.25) * (0.75 - xi - eta) * (1 - xi - eta) +
               64 * eta * (0.75 - xi - eta) * (1 - xi - eta) -
               64 * eta * (eta - 0.25) * (1 - xi - eta) -
               64 * eta * (eta - 0.25) * (0.75 - xi - eta);
    case 8:
        return 128.0 / 3 * (0.5 - xi - eta) * (0.75 - xi - eta) *
                   (1 - xi - eta) -
               128.0 / 3 * eta * (0.75 - xi - eta) * (1 - xi - eta) -
               128.0 / 3 * eta * (0.5 - xi - eta) * (1 - xi - eta) -
               128.0 / 3 * eta * (0.5 - xi - eta) * (0.75 - xi - eta);
    case 14:
        return 128 * xi * (0.75 - xi - eta) * (1 - xi - eta) -
               128 * xi * eta * (1 - xi - eta) -
               128 * xi * eta * (0.75 - xi - eta);
    case 12:
        return 128 * xi * (xi - 0.25) * (1 - xi - eta) -
               128 * xi * eta * (xi - 0.25);
    case 13:
        return 128 * xi * (eta - 0.25) * (1 - xi - eta) +
               128 * xi * eta * (1 - xi - eta) - 128 * xi * eta * (eta - 0.25);
    }
    return 0;
}

double Triangular_4_Element::det_Jacobi(double xi, double eta) {
    double J11 = 0, J12 = 0, J21 = 0, J22 = 0;
    for (int i = 0; i < n_Dofs(); i++) {
        J11 = J11 + dofslist[i][0] * phi_xi(xi, eta, i);
        J12 = J12 + dofslist[i][1] * phi_xi(xi, eta, i);
        J21 = J21 + dofslist[i][0] * phi_eta(xi, eta, i);
        J22 = J22 + dofslist[i][1] * phi_eta(xi, eta, i);
    }
    return (J11 * J22 - J12 * J21);
}

double Triangular_4_Element::xi_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++)
        a = a + dofslist[i][1] * phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_4_Element::xi_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++)
        a = a + dofslist[i][0] * phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_4_Element::eta_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++)
        a = a + dofslist[i][1] * phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_4_Element::eta_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++)
        a = a + dofslist[i][0] * phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_4_Element::phi_x(double xi, double eta, int i) {
    return phi_xi(xi, eta, i) * xi_x(xi, eta) +
           phi_eta(xi, eta, i) * eta_x(xi, eta);
}

double Triangular_4_Element::phi_y(double xi, double eta, int i) {
    return phi_xi(xi, eta, i) * xi_y(xi, eta) +
           phi_eta(xi, eta, i) * eta_y(xi, eta);
}
std::vector<double> Triangular_4_Element::gradient(double xi, double eta,
                                                   int i) {
    std::vector<double> gradient({phi_x(xi, eta, i), phi_y(xi, eta, i)});
    return gradient;
}
double Triangular_4_Element::Global_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++) {
        a += dofslist[i][0] * phi(xi, eta, i);
    }
    return a;
}

double Triangular_4_Element::Global_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < n_Dofs(); i++) {
        a += dofslist[i][1] * phi(xi, eta, i);
    }
    return a;
}

void Triangular_4_Element::reinit(Grid<2> &_grid, Mesh<2> &_mesh) {
    // TODO:if structed mesh
}
void Triangular_4_Element::reinit(Grid<2> &_grid, Mesh<2> &_mesh, int index) {
    // TODO: if unstructed mesh, you must use this function,not prevous one.
    dofslist.clear();
    int num = this->n_Dofs();
    vector<Point<2>> p = _grid.getVertices();
    vector<Point<2>> v = p;
    v.resize(15);
    // for (int j = 1; j <= 3; j++){
    //   v[3*j] = (p[j-1] * 3.0 + p[j%3])/4.0;
    //   v[3*j+1] = (p[j-1] + p[j%3])/2.0;
    //   v[3*j+2] = (p[j-1] + p[j%3] * 3.0)/4.0;
    // }
    // for(int i = 0; i < 3; i++){
    //   v[12+i] = (p[i] * 2.0 + p[(i+1)%3] + p[(i+2)%3])/4.0;
    // }
    v[3] = (v[1] * 3.0 + v[2]) / 4.0;
    v[4] = (v[1] + v[2]) / 2.0;
    v[5] = (v[1] + v[2] * 3.0) / 4.0;
    v[6] = (v[2] * 3.0 + v[0]) / 4.0;
    v[7] = (v[2] + v[0]) / 2.0;
    v[8] = (v[2] + v[0] * 3.0) / 4.0;
    v[9] = (v[0] * 3.0 + v[1]) / 4.0;
    v[10] = (v[0] + v[1]) / 2.0;
    v[11] = (v[0] + v[1] * 3.0) / 4.0;
    v[12] = (v[0] + v[1] * 2.0 + v[2]) / 4.0;
    v[13] = (v[0] + v[1] + v[2] * 2.0) / 4.0;
    v[14] = (v[0] * 2.0 + v[1] + v[2]) / 4.0;
    vector<pair<int, int>> dofsIndex = _mesh.getIndexofGrid(index);
    for (int k = 0; k < num; k++) {
        Dofs<2> q = {v[k], dofsIndex[k].first};
        dofslist.push_back(q);
    }
}

void Triangular_4_Element::reinit(Grid<2> &_grid, EasyMesh &easy_mesh,
                                  int index) {
    dofslist.clear();
    int num = this->n_Dofs();
    vector<pair<int, int>> dofsIndex = easy_mesh.getIndexofGrid(index);
    for (int k = 0; k < num; k++) {
        int index = dofsIndex[k].first;
        Dofs<2> dof = easy_mesh.getDofs(index);
        dofslist.push_back(dof);
    }
};
/**
 *@brief 一阶四边形模板单元
 */
class Quadrilateral_1_Element : public TemplateElement<2> {
    typedef vector<Dofs<2>> DofsList;
    typedef int DofsNum;

  public:
    Quadrilateral_1_Element();
    ~Quadrilateral_1_Element() = default;
    virtual Dofs<2> &getDofsList(int i);
    virtual int getGlobalIndex(int i);
    virtual Point<2> &GaussionPoint(int i) {
        return GaussInfo.GaussianPointInfo(i);
    };
    virtual double GaussionWeight(int i) {
        return GaussInfo.GaussianWeightInfo(i);
    };
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i);
    virtual double phi_eta(double xi, double eta, int i);
    virtual double det_Jacobi(double xi, double eta);
    virtual double xi_x(double xi, double eta);
    virtual double xi_y(double xi, double eta);
    virtual double eta_x(double xi, double eta);
    virtual double eta_y(double xi, double eta);
    virtual double phi_x(double xi, double eta, int i);
    virtual double phi_y(double xi, double eta, int i);
    virtual vector<double> gradient(double xi, double eta, int i);
    virtual double Global_x(double xi, double eta);
    virtual double Global_y(double xi, double eta);
    virtual double volume() { return 4 * det_Jacobi(0, 0); };
    virtual int n_Dofs() { return 4; };
    virtual int n_GaussPnt() { return GaussInfo.GaussianPointNumInfo(); };
    virtual void reinit(Grid<2> &_grid, Mesh<2> &_mesh);
    virtual void reinit(Grid<2> &_grid, Mesh<2> &_mesh, int index);
    // virtual void setMesh(Mesh<2>& _mesh);
    // virtual vector<Dofs<2>> &getBoundaryDofsonGrid();

  private:
    DofsList dofslist;         /*< 存储绑定的Grid中自由度信息*/
    DofsList boundarydofslist; /*< 存储绑定的Grid中边界自由度信息*/
    Quadrilateral_Element_GaussianInfo
        GaussInfo; /*<Gauss积分所需：积分点和积分权重*/
};

// void Quadrilateral_1_Element::setMesh(Mesh<2>& _mesh)
// {
//     mesh.getSegment() = _mesh.getSegment();
//     mesh.getGrid() = _mesh.getGrid();
// }

Quadrilateral_1_Element::Quadrilateral_1_Element() {
    dofslist.resize(4);
    GaussInfo.LoadGaussianInfo(3); // two dimensions,tensor-product rules
}

Dofs<2> &Quadrilateral_1_Element::getDofsList(int i) { return dofslist[i]; }

double Quadrilateral_1_Element::phi(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return (xi - 1) * (eta - 1) / 4.0;
    case 1:
        return -(xi + 1) * (eta - 1) / 4.0;
    case 2:
        return (xi + 1) * (eta + 1) / 4.0;
    case 3:
        return -(xi - 1) * (eta + 1) / 4.0;
    }
    return 0;
}

double Quadrilateral_1_Element::phi_xi(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return (eta - 1) / 4.0;
    case 1:
        return -(eta - 1) / 4.0;
    case 2:
        return (eta + 1) / 4.0;
    case 3:
        return -(eta + 1) / 4.0;
    }
    return 0;
}

double Quadrilateral_1_Element::phi_eta(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return (xi - 1) / 4.0;
    case 1:
        return -(xi + 1) / 4.0;
    case 2:
        return (xi + 1) / 4.0;
    case 3:
        return -(xi - 1) / 4.0;
    }
    return 0;
}

double Quadrilateral_1_Element::det_Jacobi(double xi, double eta) {
    double J11 = 0;
    double J12 = 0;
    double J21 = 0;
    double J22 = 0;
    for (int i = 0; i < 4; i++) {
        J11 = J11 + dofslist[i][0] * phi_xi(xi, eta, i);
        J12 = J12 + dofslist[i][1] * phi_xi(xi, eta, i);
        J21 = J21 + dofslist[i][0] * phi_eta(xi, eta, i);
        J22 = J22 + dofslist[i][1] * phi_eta(xi, eta, i);
    }
    return (J11 * J22 - J12 * J21);
}

double Quadrilateral_1_Element::xi_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 4; i++)
        a = a + dofslist[i][1] * phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_1_Element::xi_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 4; i++)
        a = a + dofslist[i][0] * phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_1_Element::eta_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 4; i++)
        a = a + dofslist[i][1] * phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_1_Element::eta_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 4; i++)
        a = a + dofslist[i][0] * phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_1_Element::phi_x(double xi, double eta, int i) {
    return phi_xi(xi, eta, i) * xi_x(xi, eta) +
           phi_eta(xi, eta, i) * eta_x(xi, eta);
}

double Quadrilateral_1_Element::phi_y(double xi, double eta, int i) {
    return phi_xi(xi, eta, i) * xi_y(xi, eta) +
           phi_eta(xi, eta, i) * eta_y(xi, eta);
}

vector<double> Quadrilateral_1_Element::gradient(double xi, double eta, int i) {
    vector<double> gradient({phi_x(xi, eta, i), phi_y(xi, eta, i)});
    return gradient;
}
double Quadrilateral_1_Element::Global_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 4; i++)
        a += dofslist[i][0] * phi(xi, eta, i);
    return a;
}

double Quadrilateral_1_Element::Global_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 4; i++)
        a += dofslist[i][1] * phi(xi, eta, i);
    return a;
}

int Quadrilateral_1_Element::getGlobalIndex(int i) {
    return dofslist[i].getGlobalIndex();
}

void Quadrilateral_1_Element::reinit(Grid<2> &_grid, Mesh<2> &_mesh) {
    // gridIndex should be computed before this function.
    dofslist.clear();
    vector<Point<2>> p;
    p.resize(4);
    for (int i = 0; i < 4; i++)
        p[i] = _grid.getVertices(i);
    int j = (int)round((p[0][0] + 1) * 0.5 * _mesh.getSegment()[0]);
    int k = (int)round((p[0][1] + 1) * 0.5 * _mesh.getSegment()[1]);
    int gridIndex = j + k * _mesh.getSegment()[0];
    // Point<2> p1 = {Global_x(-1, -1), Global_y(-1, -1)};
    // Point<2> p2 = {Global_x(1, -1), Global_y(1, -1)};
    // Point<2> p3 = {Global_x(1, 1), Global_y(1, 1)};
    // Point<2> p4 = {Global_x(-1, 1), Global_y(-1, 1)};
    vector<pair<int, int>> dofsList = _mesh.getIndexofGrid(gridIndex);
    for (int k = 0; k < 4; k++) {
        Dofs<2> q = {p[k], dofsList[k].first};
        dofslist.push_back(q);
        // if (_mesh.getIndexofGrid(i)[k].second != 0)
        //     boundarydofslist.push_back(q);
    }
};

void Quadrilateral_1_Element::reinit(Grid<2> &_grid, Mesh<2> &_mesh,
                                     int index) {
    // TODO: if unstructed mesh, you must use this function,not prevous one.
    dofslist.clear();
    int num = this->n_Dofs();
    vector<Point<2>> p(num);
    for (int i = 0; i < num; i++) {
        p[i] = _grid.getVertices(i);
    }
    vector<pair<int, int>> dofsIndex = _mesh.getIndexofGrid(index);
    for (int k = 0; k < num; k++) {
        Dofs<2> q = {p[k], dofsIndex[k].first};
        dofslist.push_back(q);
    }
}
// vector<Dofs<2>> &Quadrilateral_1_Element::getBoundaryDofsonGrid()
// {
//     return boundarydofslist;
// }

class Quadrilateral_2_Element : public TemplateElement<2> {
    typedef vector<Dofs<2>> DofsList;
    typedef int DofsNum;

  public:
    Quadrilateral_2_Element();
    ~Quadrilateral_2_Element(){};
    virtual Dofs<2> &getDofsList(int i);
    virtual int getGlobalIndex(int i);
    virtual Point<2> &GaussionPoint(int i) {
        return GaussInfo.GaussianPointInfo(i);
    };
    virtual double GaussionWeight(int i) {
        return GaussInfo.GaussianWeightInfo(i);
    };
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i);
    virtual double phi_eta(double xi, double eta, int i);
    virtual double det_Jacobi(double xi, double eta);
    virtual double xi_x(double xi, double eta);
    virtual double xi_y(double xi, double eta);
    virtual double eta_x(double xi, double eta);
    virtual double eta_y(double xi, double eta);
    virtual double phi_x(double xi, double eta, int i);
    virtual double phi_y(double xi, double eta, int i);
    virtual vector<double> gradient(double xi, double eta, int i);
    virtual double volume() { return 4 * det_Jacobi(0, 0); };
    virtual double Global_x(double xi, double eta);
    virtual double Global_y(double xi, double eta);
    virtual int n_Dofs() { return 9; };
    virtual int n_GaussPnt() { return GaussInfo.GaussianPointNumInfo(); };
    virtual void reinit(Grid<2> &_grid, Mesh<2> &_mesh);
    // virtual vector<Dofs<2>> &getBoundaryDofsonGrid();

  private:
    DofsList dofslist; /*< 存储绑定的Grid中自由度信息*/
    // DofsList boundarydofslist;                    /*<
    // 存储绑定的Grid中边界自由度信息*/
    Quadrilateral_Element_GaussianInfo
        GaussInfo; /*<Gauss积分所需：积分点和积分权重*/
};

Quadrilateral_2_Element::Quadrilateral_2_Element() {
    dofslist.resize(9);
    GaussInfo.LoadGaussianInfo(5);
}

Dofs<2> &Quadrilateral_2_Element::getDofsList(int i) { return dofslist[i]; }

double Quadrilateral_2_Element::phi(double xi, double eta, int i) {
    /// 注意基函数顺序与打点分布的一致性
    switch (i) {
    case 0:
        return (1 / 4.0) * xi * eta * (1.0 - xi) * (1.0 - eta);
    case 1:
        return -(1 / 2.0) * (1.0 - eta) * eta * (1.0 - xi * xi);
    case 2:
        return -(1 / 4.0) * xi * eta * (1.0 + xi) * (1.0 - eta);
    case 3:
        return (1 / 2.0) * (1.0 + xi) * xi * (1.0 - eta * eta);
    case 4:
        return (1 / 4.0) * xi * eta * (1.0 + xi) * (1.0 + eta);
    case 5:
        return (1 / 2.0) * (1.0 + eta) * eta * (1.0 - xi * xi);
    case 6:
        return -(1 / 4.0) * xi * eta * (1.0 - xi) * (1.0 + eta);
    case 7:
        return -(1 / 2.0) * (1.0 - xi) * xi * (1.0 - eta * eta);
    case 8:
        return (1.0 - xi * xi) * (1.0 - eta * eta);
    }
    return 0;
}

double Quadrilateral_2_Element::phi_xi(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return (1 / 4.0) * (1.0 - 2.0 * xi) * eta * (1.0 - eta);
    case 1:
        return xi * (1.0 - eta) * eta;
    case 2:
        return -(1 / 4.0) * (1.0 + 2.0 * xi) * eta * (1.0 - eta);
    case 3:
        return (1 / 2.0) * (1.0 + 2.0 * xi) * (1.0 - eta * eta);
    case 4:
        return (1 / 4.0) * (1.0 + 2.0 * xi) * eta * (1.0 + eta);
    case 5:
        return -xi * (1.0 + eta) * eta;
    case 6:
        return -(1 / 4.0) * (1.0 - 2.0 * xi) * eta * (1.0 + eta);
    case 7:
        return -(1 / 2.0) * (1.0 - 2.0 * xi) * (1.0 - eta * eta);
    case 8:
        return -2.0 * xi * (1.0 - eta * eta);
    }
    return 0;
}

double Quadrilateral_2_Element::phi_eta(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return (1 / 4.0) * xi * (1.0 - xi) * (1.0 - 2.0 * eta);
    case 1:
        return -(1 / 2.0) * (1.0 - xi * xi) * (1.0 - 2.0 * eta);
    case 2:
        return -(1 / 4.0) * xi * (1.0 + xi) * (1 - 2.0 * eta);
    case 3:
        return -xi * (1.0 + xi) * eta;
    case 4:
        return (1 / 4.0) * xi * (1.0 + xi) * (1.0 + 2.0 * eta);
    case 5:
        return (1 / 2.0) * (1.0 - xi * xi) * (1.0 + 2.0 * eta);
    case 6:
        return -(1 / 4.0) * xi * (1.0 - xi) * (1 + 2.0 * eta);
    case 7:
        return xi * (1.0 - xi) * eta;
    case 8:
        return -2.0 * (1.0 - xi * xi) * eta;
    }
    return 0;
}

double Quadrilateral_2_Element::xi_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 9; i++)
        a = a + dofslist[i][1] * phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_2_Element::xi_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 9; i++)
        a = a + dofslist[i][0] * phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}
double Quadrilateral_2_Element::eta_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 9; i++)
        a = a + dofslist[i][1] * phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_2_Element::eta_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 9; i++)
        a = a + dofslist[i][0] * phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_2_Element::phi_x(double xi, double eta, int i) {
    return phi_xi(xi, eta, i) * xi_x(xi, eta) +
           phi_eta(xi, eta, i) * eta_x(xi, eta);
}

double Quadrilateral_2_Element::phi_y(double xi, double eta, int i) {
    return phi_xi(xi, eta, i) * xi_y(xi, eta) +
           phi_eta(xi, eta, i) * eta_y(xi, eta);
}

vector<double> Quadrilateral_2_Element::gradient(double xi, double eta, int i) {
    vector<double> gradient({phi_x(xi, eta, i), phi_y(xi, eta, i)});
    return gradient;
}

double Quadrilateral_2_Element::det_Jacobi(double xi, double eta) {
    double J11 = 0;
    double J12 = 0;
    double J21 = 0;
    double J22 = 0;
    for (int i = 0; i < 9; i++) {
        J11 = J11 + dofslist[i][0] * phi_xi(xi, eta, i);
        J12 = J12 + dofslist[i][1] * phi_xi(xi, eta, i);
        J21 = J21 + dofslist[i][0] * phi_eta(xi, eta, i);
        J22 = J22 + dofslist[i][1] * phi_eta(xi, eta, i);
    }
    return (J11 * J22 - J12 * J21);
}

double Quadrilateral_2_Element::Global_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 9; i++)
        a += dofslist[i][0] * phi(xi, eta, i);
    // a += dofslist[0][0] * 0.25 * (1 - xi) * (1 - eta);
    // a += dofslist[2][0] * 0.25 * (1 + xi) * (1 - eta);
    // a += dofslist[4][0] * 0.25 * (1 + xi) * (1 + eta);
    // a += dofslist[6][0] * 0.25 * (1 - xi) * (1 + eta);
    return a;
}

double Quadrilateral_2_Element::Global_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 9; i++)
        a += dofslist[i][1] * phi(xi, eta, i);
    // a += dofslist[0][1] * 0.25 * (1 - xi) * (1 - eta);
    // a += dofslist[2][1] * 0.25 * (1 + xi) * (1 - eta);
    // a += dofslist[4][1] * 0.25 * (1 + xi) * (1 + eta);
    // a += dofslist[6][1] * 0.25 * (1 - xi) * (1 + eta);
    return a;
}

int Quadrilateral_2_Element::getGlobalIndex(int i) {
    return dofslist[i].getGlobalIndex();
}

void Quadrilateral_2_Element::reinit(Grid<2> &_grid, Mesh<2> &_mesh) {
    dofslist.clear();
    vector<Point<2>> p;
    p.resize(4);
    for (int i = 0; i < 4; i++)
        p[i] = _grid.getVertices(i);
    int j = (int)round((p[0][0] + 1) * 0.5 * _mesh.getSegment()[0]);
    int k = (int)round((p[0][1] + 1) * 0.5 * _mesh.getSegment()[1]);
    int gridIndex = j + k * _mesh.getSegment()[0];
    vector<Dofs<2>> q;
    q.resize(9);
    vector<pair<int, int>> dofsList = _mesh.getIndexofGrid(gridIndex);
    q[0] = {p[0], dofsList[0].first};
    q[1] = {(p[0] + p[1]) / 2, dofsList[1].first};
    q[2] = {p[1], dofsList[2].first};
    q[3] = {(p[1] + p[2]) / 2, dofsList[3].first};
    q[4] = {p[2], dofsList[4].first};
    q[5] = {(p[2] + p[3]) / 2, dofsList[5].first};
    q[6] = {p[3], dofsList[6].first};
    q[7] = {(p[3] + p[0]) / 2, dofsList[7].first};
    q[8] = {(p[0] + p[1] + p[2] + p[3]) / 4, dofsList[8].first};
    for (int k = 0; k < 9; k++) {
        dofslist.push_back(q[k]);
        // if (_mesh.getIndexofGrid(i)[k].second != 0)
        //     boundarydofslist.push_back(q[k]);
    }
}

// vector<Dofs<2>> &Quadrilateral_2_Element::getBoundaryDofsonGrid()
// {
//     return boundarydofslist;
// }

class Quadrilateral_4_Element : public TemplateElement<2> {
    typedef vector<Dofs<2>> DofsList;
    typedef int DofsNum;

  public:
    Quadrilateral_4_Element();
    ~Quadrilateral_4_Element(){};
    virtual Dofs<2> &getDofsList(int i);
    virtual int getGlobalIndex(int i);
    virtual Point<2> &GaussionPoint(int i) {
        return GaussInfo.GaussianPointInfo(i);
    };
    virtual double GaussionWeight(int i) {
        return GaussInfo.GaussianWeightInfo(i);
    };
    virtual double phi_base(double x, int i);
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i);
    virtual double phi_eta(double xi, double eta, int i);
    virtual double det_Jacobi(double xi, double eta);
    virtual double xi_x(double xi, double eta);
    virtual double xi_y(double xi, double eta);
    virtual double eta_x(double xi, double eta);
    virtual double eta_y(double xi, double eta);
    virtual double phi_x(double xi, double eta, int i);
    virtual double phi_y(double xi, double eta, int i);
    virtual vector<double> gradient(double xi, double eta, int i);
    virtual double volume() { return 4 * det_Jacobi(0, 0); };
    virtual double Global_x(double xi, double eta);
    virtual double Global_y(double xi, double eta);
    virtual int n_Dofs() { return 25; };
    virtual int n_GaussPnt() { return GaussInfo.GaussianPointNumInfo(); };
    virtual void reinit(Grid<2> &_grid, Mesh<2> &_mesh);
    // virtual vector<Dofs<2>> &getBoundaryDofsonGrid();

  private:
    DofsList dofslist; /*< 存储绑定的Grid中自由度信息*/
    Quadrilateral_Element_GaussianInfo
        GaussInfo; /*<Gauss积分所需：积分点和积分权重*/
};

Quadrilateral_4_Element::Quadrilateral_4_Element() {
    dofslist.resize(25);
    GaussInfo.LoadGaussianInfo(9);
}

Dofs<2> &Quadrilateral_4_Element::getDofsList(int i) { return dofslist[i]; }

double Quadrilateral_4_Element::phi_base(double x, int i) {
    switch (i) {
    case 0:
        return (2 / 3.0) * (x + 0.5) * x * (x - 0.5) * (x - 1);
    case 1:
        return -(8 / 3.0) * (x + 1) * x * (x - 0.5) * (x - 1);
    case 2:
        return 4 * (x + 1) * (x + 0.5) * (x - 0.5) * (x - 1);
    case 3:
        return -(8 / 3.0) * (x + 1) * (x + 0.5) * x * (x - 1);
    case 4:
        return (2 / 3.0) * (x + 1) * (x + 0.5) * x * (x - 0.5);
    }
    return 0;
}
double Quadrilateral_4_Element::phi(double xi, double eta, int i) {
    /// 注意基函数顺序与打点分布的一致性
    switch (i) {
    case 0:
        return phi_base(xi, 0) * phi_base(eta, 0);
    case 1:
        return phi_base(xi, 1) * phi_base(eta, 0);
    case 2:
        return phi_base(xi, 2) * phi_base(eta, 0);
    case 3:
        return phi_base(xi, 3) * phi_base(eta, 0);
    case 4:
        return phi_base(xi, 4) * phi_base(eta, 0);

    case 5:
        return phi_base(xi, 0) * phi_base(eta, 1);
    case 6:
        return phi_base(xi, 1) * phi_base(eta, 1);
    case 7:
        return phi_base(xi, 2) * phi_base(eta, 1);
    case 8:
        return phi_base(xi, 3) * phi_base(eta, 1);
    case 9:
        return phi_base(xi, 4) * phi_base(eta, 1);

    case 10:
        return phi_base(xi, 0) * phi_base(eta, 2);
    case 11:
        return phi_base(xi, 1) * phi_base(eta, 2);
    case 12:
        return phi_base(xi, 2) * phi_base(eta, 2);
    case 13:
        return phi_base(xi, 3) * phi_base(eta, 2);
    case 14:
        return phi_base(xi, 4) * phi_base(eta, 2);

    case 15:
        return phi_base(xi, 0) * phi_base(eta, 3);
    case 16:
        return phi_base(xi, 1) * phi_base(eta, 3);
    case 17:
        return phi_base(xi, 2) * phi_base(eta, 3);
    case 18:
        return phi_base(xi, 3) * phi_base(eta, 3);
    case 19:
        return phi_base(xi, 4) * phi_base(eta, 3);

    case 20:
        return phi_base(xi, 0) * phi_base(eta, 4);
    case 21:
        return phi_base(xi, 1) * phi_base(eta, 4);
    case 22:
        return phi_base(xi, 2) * phi_base(eta, 4);
    case 23:
        return phi_base(xi, 3) * phi_base(eta, 4);
    case 24:
        return phi_base(xi, 4) * phi_base(eta, 4);
    }
    return 0;
}

double Quadrilateral_4_Element::phi_xi(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return (4 / 9.0) * (xi * eta * (xi - 1) * (xi - 0.5) * (eta - 1) *
                                (eta - 0.5) * (eta + 0.5) +
                            xi * eta * (xi - 1) * (xi + 0.5) * (eta - 1) *
                                (eta - 0.5) * (eta + 0.5) +
                            xi * eta * (xi - 0.5) * (xi + 0.5) * (eta - 1) *
                                (eta - 0.5) * (eta + 0.5) +
                            eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (eta - 1) * (eta - 0.5) * (eta + 0.5));
    case 1:
        return (-16 / 9.0) * (xi * eta * (xi - 1) * (xi - 0.5) * (eta - 1) *
                                  (eta - 0.5) * (eta + 0.5) +
                              xi * eta * (xi - 1) * (xi + 1) * (eta - 1) *
                                  (eta - 0.5) * (eta + 0.5) +
                              xi * eta * (xi - 0.5) * (xi + 1) * (eta - 1) *
                                  (eta - 0.5) * (eta + 0.5) +
                              eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) * (eta + 0.5));
    case 2:
        return (8 / 3.0) * (eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (eta - 1) * (eta - 0.5) * (eta + 0.5) +
                            eta * (xi - 1) * (xi - 0.5) * (xi + 1) * (eta - 1) *
                                (eta - 0.5) * (eta + 0.5) +
                            eta * (xi - 1) * (xi + 0.5) * (xi + 1) * (eta - 1) *
                                (eta - 0.5) * (eta + 0.5) +
                            eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 1) * (eta - 0.5) * (eta + 0.5));
    case 3:
        return (-16 / 9.0) * (xi * eta * (xi - 1) * (xi + 0.5) * (eta - 1) *
                                  (eta - 0.5) * (eta + 0.5) +
                              xi * eta * (xi - 1) * (xi + 1) * (eta - 1) *
                                  (eta - 0.5) * (eta + 0.5) +
                              xi * eta * (xi + 0.5) * (xi + 1) * (eta - 1) *
                                  (eta - 0.5) * (eta + 0.5) +
                              eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) * (eta + 0.5));
    case 4:
        return (4 / 9.0) * (xi * eta * (xi - 0.5) * (xi + 0.5) * (eta - 1) *
                                (eta - 0.5) * (eta + 0.5) +
                            xi * eta * (xi - 0.5) * (xi + 1) * (eta - 1) *
                                (eta - 0.5) * (eta + 0.5) +
                            xi * eta * (xi + 0.5) * (xi + 1) * (eta - 1) *
                                (eta - 0.5) * (eta + 0.5) +
                            eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 1) * (eta - 0.5) * (eta + 0.5));

    case 5:
        return (-16 / 9.0) * (xi * eta * (xi - 1) * (xi - 0.5) * (eta - 1) *
                                  (eta - 0.5) * (eta + 1) +
                              xi * eta * (xi - 1) * (xi + 0.5) * (eta - 1) *
                                  (eta - 0.5) * (eta + 1) +
                              xi * eta * (xi - 0.5) * (xi + 0.5) * (eta - 1) *
                                  (eta - 0.5) * (eta + 1) +
                              eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (eta - 1) * (eta - 0.5) * (eta + 1));
    case 6:
        return (64 / 9.0) * (xi * eta * (xi - 1) * (xi - 0.5) * (eta - 1) *
                                 (eta - 0.5) * (eta + 1) +
                             xi * eta * (xi - 1) * (xi + 1) * (eta - 1) *
                                 (eta - 0.5) * (eta + 1) +
                             xi * eta * (xi - 0.5) * (xi + 1) * (eta - 1) *
                                 (eta - 0.5) * (eta + 1) +
                             eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                 (eta - 1) * (eta - 0.5) * (eta + 1));
    case 7:
        return (-32 / 3.0) * (eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (eta - 1) * (eta - 0.5) * (eta + 1) +
                              eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) * (eta + 1) +
                              eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) * (eta + 1) +
                              eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) * (eta + 1));
    case 8:
        return (64 / 9.0) * (xi * eta * (xi - 1) * (xi + 0.5) * (eta - 1) *
                                 (eta - 0.5) * (eta + 1) +
                             xi * eta * (xi - 1) * (xi + 1) * (eta - 1) *
                                 (eta - 0.5) * (eta + 1) +
                             xi * eta * (xi + 0.5) * (xi + 1) * (eta - 1) *
                                 (eta - 0.5) * (eta + 1) +
                             eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                 (eta - 1) * (eta - 0.5) * (eta + 1));
    case 9:
        return (-16 / 9.0) * (xi * eta * (xi - 0.5) * (xi + 0.5) * (eta - 1) *
                                  (eta - 0.5) * (eta + 1) +
                              xi * eta * (xi - 0.5) * (xi + 1) * (eta - 1) *
                                  (eta - 0.5) * (eta + 1) +
                              xi * eta * (xi + 0.5) * (xi + 1) * (eta - 1) *
                                  (eta - 0.5) * (eta + 1) +
                              eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) * (eta + 1));

    case 10:
        return (8 / 3.0) * (xi * (xi - 1) * (xi - 0.5) * (eta - 1) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                            xi * (xi - 1) * (xi + 0.5) * (eta - 1) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                            xi * (xi - 0.5) * (xi + 0.5) * (eta - 1) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                            (xi - 1) * (xi - 0.5) * (xi + 0.5) * (eta - 1) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1));
    case 11:
        return (-32 / 3.0) * (xi * (xi - 1) * (xi - 0.5) * (eta - 1) *
                                  (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                              xi * (xi - 1) * (xi + 1) * (eta - 1) *
                                  (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                              xi * (xi - 0.5) * (xi + 1) * (eta - 1) *
                                  (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                              (xi - 1) * (xi - 0.5) * (xi + 1) * (eta - 1) *
                                  (eta - 0.5) * (eta + 0.5) * (eta + 1));
    case 12:
        return 16 * ((xi - 1) * (xi - 0.5) * (xi + 0.5) * (eta - 1) *
                         (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                     (xi - 1) * (xi - 0.5) * (xi + 1) * (eta - 1) *
                         (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                     (xi - 1) * (xi + 0.5) * (xi + 1) * (eta - 1) *
                         (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                     (xi - 0.5) * (xi + 0.5) * (xi + 1) * (eta - 1) *
                         (eta - 0.5) * (eta + 0.5) * (eta + 1));
    case 13:
        return (-32 / 3.0) * (xi * (xi - 1) * (xi + 0.5) * (eta - 1) *
                                  (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                              xi * (xi - 1) * (xi + 1) * (eta - 1) *
                                  (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                              xi * (xi + 0.5) * (xi + 1) * (eta - 1) *
                                  (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                              (xi - 1) * (xi + 0.5) * (xi + 1) * (eta - 1) *
                                  (eta - 0.5) * (eta + 0.5) * (eta + 1));
    case 14:
        return (8 / 3.0) * (xi * (xi - 0.5) * (xi + 0.5) * (eta - 1) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                            xi * (xi - 0.5) * (xi + 1) * (eta - 1) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                            xi * (xi + 0.5) * (xi + 1) * (eta - 1) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                            (xi - 0.5) * (xi + 0.5) * (xi + 1) * (eta - 1) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1));

    case 15:
        return (-16 / 9.0) * (xi * eta * (xi - 1) * (xi - 0.5) * (eta - 1) *
                                  (eta + 0.5) * (eta + 1) +
                              xi * eta * (xi - 1) * (xi + 0.5) * (eta - 1) *
                                  (eta + 0.5) * (eta + 1) +
                              xi * eta * (xi - 0.5) * (xi + 0.5) * (eta - 1) *
                                  (eta + 0.5) * (eta + 1) +
                              eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (eta - 1) * (eta + 0.5) * (eta + 1));
    case 16:
        return (64 / 9.0) * (xi * eta * (xi - 1) * (xi - 0.5) * (eta - 1) *
                                 (eta + 0.5) * (eta + 1) +
                             xi * eta * (xi - 1) * (xi + 1) * (eta - 1) *
                                 (eta + 0.5) * (eta + 1) +
                             xi * eta * (xi - 0.5) * (xi + 1) * (eta - 1) *
                                 (eta + 0.5) * (eta + 1) +
                             eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                 (eta - 1) * (eta + 0.5) * (eta + 1));
    case 17:
        return (-32 / 3.0) * (eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (eta - 1) * (eta + 0.5) * (eta + 1) +
                              eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                  (eta - 1) * (eta + 0.5) * (eta + 1) +
                              eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta + 0.5) * (eta + 1) +
                              eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta + 0.5) * (eta + 1));
    case 18:
        return (64 / 9.0) * (xi * eta * (xi - 1) * (xi + 0.5) * (eta - 1) *
                                 (eta + 0.5) * (eta + 1) +
                             xi * eta * (xi - 1) * (xi + 1) * (eta - 1) *
                                 (eta + 0.5) * (eta + 1) +
                             xi * eta * (xi + 0.5) * (xi + 1) * (eta - 1) *
                                 (eta + 0.5) * (eta + 1) +
                             eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                 (eta - 1) * (eta + 0.5) * (eta + 1));
    case 19:
        return (-16 / 9.0) * (xi * eta * (xi - 0.5) * (xi + 0.5) * (eta - 1) *
                                  (eta + 0.5) * (eta + 1) +
                              xi * eta * (xi - 0.5) * (xi + 1) * (eta - 1) *
                                  (eta + 0.5) * (eta + 1) +
                              xi * eta * (xi + 0.5) * (xi + 1) * (eta - 1) *
                                  (eta + 0.5) * (eta + 1) +
                              eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta + 0.5) * (eta + 1));

    case 20:
        return (4 / 9.0) * (xi * eta * (xi - 1) * (xi - 0.5) * (eta - 0.5) *
                                (eta + 0.5) * (eta + 1) +
                            xi * eta * (xi - 1) * (xi + 0.5) * (eta - 0.5) *
                                (eta + 0.5) * (eta + 1) +
                            xi * eta * (xi - 0.5) * (xi + 0.5) * (eta - 0.5) *
                                (eta + 0.5) * (eta + 1) +
                            eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1));
    case 21:
        return (-16 / 9.0) * (xi * eta * (xi - 1) * (xi - 0.5) * (eta - 0.5) *
                                  (eta + 0.5) * (eta + 1) +
                              xi * eta * (xi - 1) * (xi + 1) * (eta - 0.5) *
                                  (eta + 0.5) * (eta + 1) +
                              xi * eta * (xi - 0.5) * (xi + 1) * (eta - 0.5) *
                                  (eta + 0.5) * (eta + 1) +
                              eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                  (eta - 0.5) * (eta + 0.5) * (eta + 1));
    case 22:
        return (8 / 3.0) * (eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                            eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                            eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1) +
                            eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1));
    case 23:
        return (-16 / 9.0) * (xi * eta * (xi - 1) * (xi + 0.5) * (eta - 0.5) *
                                  (eta + 0.5) * (eta + 1) +
                              xi * eta * (xi - 1) * (xi + 1) * (eta - 0.5) *
                                  (eta + 0.5) * (eta + 1) +
                              xi * eta * (xi + 0.5) * (xi + 1) * (eta - 0.5) *
                                  (eta + 0.5) * (eta + 1) +
                              eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                  (eta - 0.5) * (eta + 0.5) * (eta + 1));
    case 24:
        return (4 / 9.0) * (xi * eta * (xi - 0.5) * (xi + 0.5) * (eta - 0.5) *
                                (eta + 0.5) * (eta + 1) +
                            xi * eta * (xi - 0.5) * (xi + 1) * (eta - 0.5) *
                                (eta + 0.5) * (eta + 1) +
                            xi * eta * (xi + 0.5) * (xi + 1) * (eta - 0.5) *
                                (eta + 0.5) * (eta + 1) +
                            eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1));
    }
    return 0;
}

double Quadrilateral_4_Element::phi_eta(double xi, double eta, int i) {
    switch (i) {
    case 0:
        return (4 / 9.0) * (xi * eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (eta - 1) * (eta - 0.5) +
                            xi * eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (eta - 1) * (eta + 0.5) +
                            xi * eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (eta - 0.5) * (eta + 0.5) +
                            xi * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (eta - 1) * (eta - 0.5) * (eta + 0.5));
    case 1:
        return (-16 / 9.0) * (xi * eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) +
                              xi * eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                  (eta - 1) * (eta + 0.5) +
                              xi * eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                  (eta - 0.5) * (eta + 0.5) +
                              xi * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) * (eta + 0.5));
    case 2:
        return (8 / 3.0) * (eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (xi + 1) * (eta - 1) * (eta - 0.5) +
                            eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (xi + 1) * (eta - 1) * (eta + 0.5) +
                            eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (xi + 1) * (eta - 0.5) * (eta + 0.5) +
                            (xi - 1) * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 1) * (eta - 0.5) * (eta + 0.5));
    case 3:
        return (-16 / 9.0) * (xi * eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) +
                              xi * eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta + 0.5) +
                              xi * eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                  (eta - 0.5) * (eta + 0.5) +
                              xi * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) * (eta + 0.5));
    case 4:
        return (4 / 9.0) * (xi * eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 1) * (eta - 0.5) +
                            xi * eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 1) * (eta + 0.5) +
                            xi * eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 0.5) * (eta + 0.5) +
                            xi * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 1) * (eta - 0.5) * (eta + 0.5));

    case 5:
        return (-16 / 9.0) * (xi * eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (eta - 1) * (eta - 0.5) +
                              xi * eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (eta - 1) * (eta + 1) +
                              xi * eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (eta - 0.5) * (eta + 1) +
                              xi * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (eta - 1) * (eta - 0.5) * (eta + 1));
    case 6:
        return (64 / 9.0) * (xi * eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                 (eta - 1) * (eta - 0.5) +
                             xi * eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                 (eta - 1) * (eta + 1) +
                             xi * eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                 (eta - 0.5) * (eta + 1) +
                             xi * (xi - 1) * (xi - 0.5) * (xi + 1) * (eta - 1) *
                                 (eta - 0.5) * (eta + 1));
    case 7:
        return (-32 / 3.0) * (eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (xi + 1) * (eta - 1) * (eta - 0.5) +
                              eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (xi + 1) * (eta - 1) * (eta + 1) +
                              eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (xi + 1) * (eta - 0.5) * (eta + 1) +
                              (xi - 1) * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) * (eta + 1));
    case 8:
        return (64 / 9.0) * (xi * eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                 (eta - 1) * (eta - 0.5) +
                             xi * eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                 (eta - 1) * (eta + 1) +
                             xi * eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                 (eta - 0.5) * (eta + 1) +
                             xi * (xi - 1) * (xi + 0.5) * (xi + 1) * (eta - 1) *
                                 (eta - 0.5) * (eta + 1));
    case 9:
        return (-16 / 9.0) * (xi * eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) +
                              xi * eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta + 1) +
                              xi * eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                  (eta - 0.5) * (eta + 1) +
                              xi * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) * (eta + 1));

    case 10:
        return (8 / 3.0) * (xi * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (eta - 1) * (eta - 0.5) * (eta + 0.5) +
                            xi * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (eta - 1) * (eta - 0.5) * (eta + 1) +
                            xi * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (eta - 1) * (eta + 0.5) * (eta + 1) +
                            xi * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1));
    case 11:
        return (-32 / 3.0) * (xi * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) * (eta + 0.5) +
                              xi * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) * (eta + 1) +
                              xi * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                  (eta - 1) * (eta + 0.5) * (eta + 1) +
                              xi * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                  (eta - 0.5) * (eta + 0.5) * (eta + 1));
    case 12:
        return 16 * ((xi - 1) * (xi - 0.5) * (xi + 0.5) * (xi + 1) * (eta - 1) *
                         (eta - 0.5) * (eta + 0.5) +
                     (xi - 1) * (xi - 0.5) * (xi + 0.5) * (xi + 1) * (eta - 1) *
                         (eta - 0.5) * (eta + 1) +
                     (xi - 1) * (xi - 0.5) * (xi + 0.5) * (xi + 1) * (eta - 1) *
                         (eta + 0.5) * (eta + 1) +
                     (xi - 1) * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                         (eta - 0.5) * (eta + 0.5) * (eta + 1));
    case 13:
        return (-32 / 3.0) * (xi * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) * (eta + 0.5) +
                              xi * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta - 0.5) * (eta + 1) +
                              xi * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta + 0.5) * (eta + 1) +
                              xi * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                  (eta - 0.5) * (eta + 0.5) * (eta + 1));
    case 14:
        return (8 / 3.0) * (xi * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 1) * (eta - 0.5) * (eta + 0.5) +
                            xi * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 1) * (eta - 0.5) * (eta + 1) +
                            xi * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 1) * (eta + 0.5) * (eta + 1) +
                            xi * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1));

    case 15:
        return (-16 / 9.0) * (xi * eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (eta - 1) * (eta + 0.5) +
                              xi * eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (eta - 1) * (eta + 1) +
                              xi * eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (eta + 0.5) * (eta + 1) +
                              xi * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (eta - 1) * (eta + 0.5) * (eta + 1));
    case 16:
        return (64 / 9.0) * (xi * eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                 (eta - 1) * (eta + 0.5) +
                             xi * eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                 (eta - 1) * (eta + 1) +
                             xi * eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                 (eta + 0.5) * (eta + 1) +
                             xi * (xi - 1) * (xi - 0.5) * (xi + 1) * (eta - 1) *
                                 (eta + 0.5) * (eta + 1));
    case 17:
        return (-32 / 3.0) * (eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (xi + 1) * (eta - 1) * (eta + 0.5) +
                              eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (xi + 1) * (eta - 1) * (eta + 1) +
                              eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                  (xi + 1) * (eta + 0.5) * (eta + 1) +
                              (xi - 1) * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta + 0.5) * (eta + 1));
    case 18:
        return (64 / 9.0) * (xi * eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                 (eta - 1) * (eta + 0.5) +
                             xi * eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                 (eta - 1) * (eta + 1) +
                             xi * eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                 (eta + 0.5) * (eta + 1) +
                             xi * (xi - 1) * (xi + 0.5) * (xi + 1) * (eta - 1) *
                                 (eta + 0.5) * (eta + 1));
    case 19:
        return (-16 / 9.0) * (xi * eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta + 0.5) +
                              xi * eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta + 1) +
                              xi * eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                  (eta + 0.5) * (eta + 1) +
                              xi * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                  (eta - 1) * (eta + 0.5) * (eta + 1));

    case 20:
        return (4 / 9.0) * (xi * eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (eta - 0.5) * (eta + 0.5) +
                            xi * eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (eta - 0.5) * (eta + 1) +
                            xi * eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (eta + 0.5) * (eta + 1) +
                            xi * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1));
    case 21:
        return (-16 / 9.0) * (xi * eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                  (eta - 0.5) * (eta + 0.5) +
                              xi * eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                  (eta - 0.5) * (eta + 1) +
                              xi * eta * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                  (eta + 0.5) * (eta + 1) +
                              xi * (xi - 1) * (xi - 0.5) * (xi + 1) *
                                  (eta - 0.5) * (eta + 0.5) * (eta + 1));
    case 22:
        return (8 / 3.0) * (eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (xi + 1) * (eta - 0.5) * (eta + 0.5) +
                            eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (xi + 1) * (eta - 0.5) * (eta + 1) +
                            eta * (xi - 1) * (xi - 0.5) * (xi + 0.5) *
                                (xi + 1) * (eta + 0.5) * (eta + 1) +
                            (xi - 1) * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1));
    case 23:
        return (-16 / 9.0) * (xi * eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                  (eta - 0.5) * (eta + 0.5) +
                              xi * eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                  (eta - 0.5) * (eta + 1) +
                              xi * eta * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                  (eta + 0.5) * (eta + 1) +
                              xi * (xi - 1) * (xi + 0.5) * (xi + 1) *
                                  (eta - 0.5) * (eta + 0.5) * (eta + 1));
    case 24:
        return (4 / 9.0) * (xi * eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 0.5) * (eta + 0.5) +
                            xi * eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 0.5) * (eta + 1) +
                            xi * eta * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta + 0.5) * (eta + 1) +
                            xi * (xi - 0.5) * (xi + 0.5) * (xi + 1) *
                                (eta - 0.5) * (eta + 0.5) * (eta + 1));
    }
    return 0;
}

double Quadrilateral_4_Element::xi_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 25; i++)
        a = a + dofslist[i][1] * phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_4_Element::xi_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 25; i++)
        a = a + dofslist[i][0] * phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}
double Quadrilateral_4_Element::eta_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 25; i++)
        a = a + dofslist[i][1] * phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_4_Element::eta_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 25; i++)
        a = a + dofslist[i][0] * phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_4_Element::phi_x(double xi, double eta, int i) {
    return phi_xi(xi, eta, i) * xi_x(xi, eta) +
           phi_eta(xi, eta, i) * eta_x(xi, eta);
}

double Quadrilateral_4_Element::phi_y(double xi, double eta, int i) {
    return phi_xi(xi, eta, i) * xi_y(xi, eta) +
           phi_eta(xi, eta, i) * eta_y(xi, eta);
}

vector<double> Quadrilateral_4_Element::gradient(double xi, double eta, int i) {
    vector<double> gradient({phi_x(xi, eta, i), phi_y(xi, eta, i)});
    return gradient;
}

double Quadrilateral_4_Element::det_Jacobi(double xi, double eta) {
    double J11 = 0;
    double J12 = 0;
    double J21 = 0;
    double J22 = 0;
    for (int i = 0; i < 25; i++) {
        J11 = J11 + dofslist[i][0] * phi_xi(xi, eta, i);
        J12 = J12 + dofslist[i][1] * phi_xi(xi, eta, i);
        J21 = J21 + dofslist[i][0] * phi_eta(xi, eta, i);
        J22 = J22 + dofslist[i][1] * phi_eta(xi, eta, i);
    }
    return (J11 * J22 - J12 * J21);
}

double Quadrilateral_4_Element::Global_x(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 25; i++)
        a += dofslist[i][0] * phi(xi, eta, i);
    // a += dofslist[0][0] * 0.25 * (1 - xi) * (1 - eta);
    // a += dofslist[4][0] * 0.25 * (1 + xi) * (1 - eta);
    // a += dofslist[24][0] * 0.25 * (1 + xi) * (1 + eta);
    // a += dofslist[20][0] * 0.25 * (1 - xi) * (1 + eta);
    return a;
}

double Quadrilateral_4_Element::Global_y(double xi, double eta) {
    double a = 0;
    for (int i = 0; i < 25; i++)
        a += dofslist[i][1] * phi(xi, eta, i);
    // a += dofslist[0][0] * 0.25 * (1 - xi) * (1 - eta);
    // a += dofslist[4][0] * 0.25 * (1 + xi) * (1 - eta);
    // a += dofslist[24][0] * 0.25 * (1 + xi) * (1 + eta);
    // a += dofslist[20][0] * 0.25 * (1 - xi) * (1 + eta);
    return a;
}

int Quadrilateral_4_Element::getGlobalIndex(int i)
{
  return dofslist[i].getGlobalIndex();
}

void Quadrilateral_4_Element::reinit(Grid<2> &_grid, Mesh<2> &_mesh)
{
  dofslist.clear();
  vector<Point<2>> p;
  p.resize(4);
  for (int i = 0; i < 4; i++)
    p[i] = _grid.getVertices(i);
  int j = (int)round((p[0][0] + 1) * 0.5 * _mesh.getSegment()[0]);
  int k = (int)round((p[0][1] + 1) * 0.5 * _mesh.getSegment()[1]);
  int gridIndex = j + k * _mesh.getSegment()[0];
  vector<Dofs<2>> q;
  q.resize(25);
  vector<pair<int, int>> dofsList = _mesh.getIndexofGrid(gridIndex);
  for (int j = 0; j < 5; j++)
  {
    for (int i = 0; i < 5; i++)
    {
      q[i + j * 5] = {(p[0] * (1 - i * 0.25) + p[1] * i * 0.25) * (1 - j * 0.25) +
        (p[3] * (1 - i * 0.25) + p[2] * i * 0.25) * j * 0.25,
          dofsList[i + j * 5].first};
      dofslist.push_back(q[i + j * 5]);
      // cout << q[i + j * 5] << endl;
    }
  }
}

#else
// Do Nothing.
#endif
